References:
#Calculating metrics
#Sources: #https://www.r-bloggers.com/how-to-calculate-landscape-metrics-for-local-landscapes/ #https://r-spatialecology.github.io/landscapemetrics/reference/list_lsm.html #https://r-spatialecology.github.io/landscapemetrics/ #https://r-spatialecology.github.io/landscapemetrics/articles/articles/guide_sample_lsm.html#extract-landscape-metrics-at-sample-points #https://r-spatialecology.github.io/landscapemetrics/reference/sample_lsm.html#details
Setting up the package with the packages we plan to use
Brings in the CSVS with the buffer as a dataframe and runs a quick test to select columns from the dataframe displaying the NDSI values per buffer per site to show the values.
Bring in the shapefile for the points of interest and examine the areas
Think about RGEE and how it could be used
The parameterized model is run on the input layers using the CAPS software, written at UMass by Brad Compton and Eduard Ene. This software produces an output grid for each metric. Metrics fall into two groups: stressor metrics (such as road traffic, invasive plants, or nutrient enrichment), and resiliency metrics (similarity, connectedness, and aquatic connectedness).
Stressor metrics measure anthropogenic stressors that reduce the integrity of a site, while resiliency metrics measure the intrinsic ability of a site to maintain its ecological integrity, despite the impact of anthropogenic stressors.
Resiliency metrics, in reflecting the currentlandscape, do take into account anthropogenic stressors such as road traffic and impervious surfaces. The three resiliency metrics are based on the ecological distance among cells computed using the ecological settings variables described in Appendix D.
Scaling. Scaled metrics and IEI are scaled from 0-1; in geoTIFFs, these grids are expressed in terms of percent (scaled 0-100). Raw metrics and settings grids are scaled in original units, unique to each grid; in geoTIFFs, these grids are scaled from 0-255. The CAPS final landcover, capsland, represents landcover classes using integer classes (see Appendix H). GeoTIFF versions are already colored appropriately; legends files for QGIS, ArcView 3.3, and ArcMap are supplied for the Arc grid.
The coordinate reference system for all data is Massachusetts mainland State Plane, NAD83.
Hydrological alterations
Imperviousness imperv Measures the intensity of impervious surface in the watershed above the focal cell, based on imperviousness and the modeled “influence value” for each cell, which is the aquatic distance from the focal cell based on a timeof-flow model. Data source: landcover, streams, flow direction, watershed resistance, percent imperviousness
Integrity Metrics
Connectedness connect Measures the disruption of habitat connectivity caused by all forms of development between each focal cell and surrounding cells as well as the “resistance” of the surrounding undeveloped landscape, as well as the similarity of surroundings. A hypothetical organism in a highly connected cell can reach a large area of ecologically similar cells with minimal crossing of “hostile” cells. This metric uses a least-cost path algorithm to determine the area that can reach each focal cell, incorporating each cell’s similarity to the focal cell. Data source: landcover, ecological settings variables
Hydrological alterations
Imperviousness
raw: http://jamba.provost.ads.umass.edu/web/CAPS2011/tiffzips/metricsraw/imperv.zip
Resiliency Metrics
Connectedness
raw:http://jamba.provost.ads.umass.edu/web/CAPS2011/tiffzips/metricsraw/connect.zip
Vegetation
Vegetative structure
http://jamba.provost.ads.umass.edu/web/caps2011/tiffzips/settings/structure.zip
Development & roads
Road traffic
raw:http://jamba.provost.ads.umass.edu/web/CAPS2011/tiffzips/metricsraw/traffic.zip
##Two ways of bringing in the Tif Files
[1] One by One
[2] All at once
###[1]
#import immperviousness raster file
imperv <- raster("../external/data/Amherst_CAPS2011/imperv.tif", crs = crs(sample_points_f))
##Plot code to test tif
# plot(imperv)
# plot(sample_points_f, col = "red", add = TRUE)
# plot(ma_crop, add = TRUE)
# st_crs(imperv)
#Crop Imperv to AOI and plot (code commented out)
imperv_cropped <- crop(imperv, y = ma_crop)
# plot(imperv_cropped)
# plot(sample_points_f, col = "red", add = TRUE)
#import connectedness raster file
connect <- raster("../external/data/Amherst_CAPS2011/connect.tif", crs = crs(sample_points_f))
##Plot code to test tif
# plot(connect)
# plot(sample_points_f, col = "red", add = TRUE)
# plot(ma_crop, add = TRUE)
# st_crs(connect)
#Crop connect to AOI and plot (code commented out)
connect_cropped <- crop(connect, y = ma_crop)
# plot(connect_cropped)
# plot(sample_points_f, col = "red", add = TRUE)
#import structure raster file
structure <- raster("../external/data/Amherst_CAPS2011/structure.tif", crs = crs(sample_points_f))
##Plot code to test tif
# plot(structure)
# plot(sample_points_f, col = "red", add = TRUE)
# plot(AOI, add = TRUE)
# st_crs(structure)
#Crop structure to AOI and plot (code commented out)
structure_cropped <- crop(structure, y = ma_crop)
# plot(structure_cropped)
# plot(sample_points_f, col = "red", add = TRUE)
# #Road traffic traffic Measures the intensity of road traffic (based on measured road traffic rates) in the neighborhood surrounding the focal cell, based on a logistic function of distance.
# Data source: landcover, traffic rates
#import connectedness raster file
traffic <- raster("../external/data/Amherst_CAPS2011/traffic.tif", crs = crs(sample_points_f))
##Plot code to test tif
# plot(traffic)
# plot(sample_points_f, col = "green", add = TRUE)
# plot(AOI, add = TRUE)
# st_crs(traffic)
#Crop connect to AOI and plot (code commented out)
traffic_cropped <- crop(traffic, y = ma_crop)
# plot(traffic_cropped)
# plot(sample_points_f, col = "green", add = TRUE)
###Example plot of data
#
# par(mar = c(0, 0, 1, 0) + .1)
# plot(structure, c, axes = FALSE, nr = 4)
# plot(sample_points_f, col = "red", add=TRUE)
#
# #Example of plotting out a single site
# par(mar = c(0, 0, 1, 0) + .1)
# plot(structure, c, axes = FALSE, nr = 4)
# sample_points_f %>% slice(1) %>% plot(col = "red", add=TRUE)
# Loading all the files together at once into a list
f <- dir("../external/data/Amherst_CAPS2011", full.names = TRUE, pattern = ".tif")
landcover_list <- list(connect_cropped, imperv_cropped, structure_cropped, traffic_cropped)
names(landcover_list) <- gsub("tif", "", basename(f))
#[2] all together
# f <- dir("./data/Amherst_CAPS2011", full.names = TRUE, pattern = ".tif")
# landcover <- lapply(unname(f), function(x) {
# r <- raster(x)
# crs(r) <- crs(pts)
# r
# })
# names(landcover) <- gsub(".tif", "", basename(f))
#Scale the images subtracting the min values of each (the zero values) and then dividing by the max value possible in the images, 255.
lc_scaled <-lapply(landcover_list,function(x){
(x- cellStats(x, min)) / 255
})
lc_stack <- stack(lc_scaled)
#Plot out all of the images at once
plot(lc_stack)Calculate focal mean stat at a weight of 3000 for each variable of interest–imperviousness, traffic, and connectedness.
#Sources for understanding focal window
# http://r-sig-geo.2731867.n2.nabble.com/Raster-package-Focal-sum-in-circles-td7587683.html
# https://gis.stackexchange.com/questions/292409/why-different-results-of-mean-calculations-with-focal-in-r-and-esri-arcgis
# https://www.timassal.com/?p=2092
# https://gis.stackexchange.com/questions/287553/how-to-use-sum-of-circular-moving-window-for-each-center-cell-with-focal-in-r
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/focal
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/focalWeight
# NOTES ON FOCAL FUNCTION
#.............................................
#focal passes over each image pixel, and multiplies those weights by each pixel value in the neighborhood, and then sums those to get the mean
# It sums the values because sum is the default value of the argument “fun” in the function focal, which is why we have not even specified the argument “fun” in Block 1
# The focal function operates on the matrix, regardless of it being circular or rectangular and assigns the resulting value to the center cell. In a binary matrix values with 1 would be operated on and those with 0, ignored. A circular window is a rectangular matrix where, values occuring outside the defined radius are 0.
# Since we are interested in the range of the recording device and because we have assigned a consistent value across the buffers for NDSI, ANT, and BIO varying per site per metric, we are going to use a focalWeight circle to replicate the 3000m buffer and pass it over the entirety of the images for imperv, connect, and structure. This will give us a raster for each and
#2 apply circular moving window to continuous data
###
#set the focal weight, since we are using a circle, set number to the radius of the circle (in units of CRS)
#which you can conventientyly check using the check_landscape() function
# #cell resolution is 30 x 30
#...................................................................
#Now let's apply the focalweight function to get a focal matrix for each .tif file in the landcover_list
#1500 Buffers
fw_sum_1500 <- lapply(landcover_list,function(x){
fw <- focalWeight(x, 1500, type='circle')
})
#2000 Buffers
fw_sum_2000 <- lapply(landcover_list,function(x){
fw <- focalWeight(x, 2000, type='circle')
})
#2500 Buffers
fw_sum_2500 <- lapply(landcover_list,function(x){
fw <- focalWeight(x, 2500, type='circle')
})
#Runs focal at a buffer weight of 1500 m radius
radius1 <- 1500
lc1500 <- lapply(1:4, function(x) { # x <- 1
f <- paste0("../external/data/Amherst_CAPS2011", names(lc_stack)[x], "_", radius1, ".tif")
r <- focal(lc_stack[[x]], w = as.matrix(fw_sum_1500[[x]]), filename = f, overwrite = TRUE)
return(r)
})
lc1500_stacked <- stack(lc1500) # stack
# Plotting out the results
tit1 <- "1500 (m) Buffer "
par(mfrow = c(2, 2), mar = c(.5, 0, 2, 6))
lapply(1:4, function(x){
plot(lc1500_stacked[[x]], axes = FALSE, box = FALSE, main = toupper(c(tit1,gsub(".tif", "", basename(f[x])))))
plot(st_geometry(sample_points_f), add = TRUE)
})## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
#Runs focal at a buffer weight of 2000 m radius
radius2 <- 2000
lc2000 <- lapply(1:4, function(x) { # x <- 1
f <- paste0("../external/data/Amherst_CAPS2011", names(lc_stack)[x], "_", radius2, ".tif")
r <- focal(lc_stack[[x]], w = as.matrix(fw_sum_2000[[x]]), filename = f, overwrite = TRUE)
return(r)
})
lc2000_stacked <- stack(lc2000) # stack
# Plotting out the results
tit2 <- "2000 (m) Buffer "
par(mfrow = c(2, 2), mar = c(.5, 0, 2, 6))
lapply(1:4, function(x){
plot(lc2000_stacked[[x]], axes = FALSE, box = FALSE, main = toupper(c(tit2,gsub(".tif", "", basename(f[x])))))
plot(st_geometry(sample_points_f), add = TRUE)
})## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
#Runs focal at a buffer weight of 2500 m radius
radius3 <- 2500
lc2500 <- lapply(1:4, function(x) { # x <- 1
f <- paste0("../external/data/Amherst_CAPS2011", names(lc_stack)[x], "_", radius3, ".tif")
r <- focal(lc_stack[[x]], w = as.matrix(fw_sum_2500[[x]]), filename = f, overwrite = TRUE)
return(r)
})
lc2500_stacked <- stack(lc2500) # stack
# Plotting out the results
tit3 <- "2500 (m) Buffer"
par(mfrow = c(2, 2), mar = c(.5, 0, 2, 6))
lapply(1:4, function(x){
plot(lc2500_stacked[[x]], axes = FALSE, box = FALSE, main = toupper(c(tit3, gsub(".tif", "", basename(f[x])))))
plot(st_geometry(sample_points_f), add = TRUE)
})## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
Write out the raster to use with the random forest prediction
# Code to write out rasters for use with the randomforest prediction
path_out <- '../inst/extdata'
writeRaster(lc2500_stacked[[1]], file.path(path_out,"lc2500_connect.tif"), overwrite = TRUE)
writeRaster(lc2500_stacked[[2]], file.path(path_out,"lc2500_imperv.tif"), overwrite = TRUE)
writeRaster(lc2500_stacked[[3]], file.path(path_out,"lc2500_structure.tif"), overwrite = TRUE)
writeRaster(lc2500_stacked[[4]], file.path(path_out, "lc2500_traffic.tif"), overwrite = TRUE)##Extract focal values to get predictors for points
pts1 <- cbind(sample_points_f, raster::extract(lc_stack, sample_points_f)) # original lc
pts_1500 <- cbind(pts1, raster::extract(lc1500_stacked, pts1))# buffer points 1500
write.csv(pts_1500, file.path(path_out,"fv_pts_1500.csv"), row.names = FALSE)
pts_2000 <- cbind(pts1, raster::extract(lc2000_stacked, pts1)) # buffer points 2000
write.csv(pts_2000, file.path(path_out,"fv_pts_2000.csv"), row.names = FALSE)
pts_2500 <- cbind(pts1, raster::extract(lc2500_stacked, pts1)) # buffer points 2500
write.csv(pts_2500, file.path(path_out,"fv_pts_2500.csv"), row.names = FALSE)Extract values from a Raster* object at the locations of other spatial data. You can use coordinates (points), lines, polygons or an Extent (rectangle) object. You can also use cell numbers to extract values. If y represents points, extract returns the values of a Raster* object for the cells in which a set of points fall. If y represents lines, the extract method returns the values of the cells of a Raster object that are touched by a line. If y represents polygons, the extract method returns the values of the cells of a Raster object that are covered by a polygon. A cell is covered if its center is inside the polygon (but see the weights option for considering partly covered cells; and argument small for getting values for small polygons). It is also possible to extract values for point locations from SpatialPolygons.
Code to create linear correlations between NDSI and IMperviousness, connectedness, and structure
Code to create linear correlations between ANT and IMperviousness, connectedness, and Structure
Code to create linear correlations between BIO and Imperviousness, connectedness, and structure
Code to create linear r squared dataframe from the rsquared summaries above
# as. Buffer <- Imperv_connect_Struct_Traffic_NDSI$Buffer
Buffer <- c(500, 1000,1500,2000,2500,3000)
linear_rsquared <- data.frame(Buffer, NDSI_Imperv_lm, NDSI_Connect_lm, NDSI_Structure_lm, NDSI_Traffic_lm, ANT_Imperv_lm, ANT_Connect_lm, ANT_Structure_lm, ANT_Traffic_lm, BIO_Imperv_lm, BIO_Connect_lm, BIO_Structure_lm, BIO_Traffic_lm)
# write.csv(linear_rsquared, file = "linear_rsquared_traffic_included.csv")
#-------------------------------------------------------------------------------------------------------------------------
#Graphs
#graph of R^2 for NDSI
RSquaredPlot_NDSI <- ggplot(data = linear_rsquared) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$NDSI_Imperv_lm, color = 'Imperv')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$NDSI_Structure_lm, color = 'Structure')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$NDSI_Connect_lm, color = 'Connect')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$NDSI_Traffic_lm, color = 'Traffic')) + scale_x_log10() +
labs(title = "Mean NDSI ~ R squared of metrics", x ="Focal Distance (meters)", y = "R squared of metrics")
plot(RSquaredPlot_NDSI)#graph of R^2 for ANT
RSquaredPlot_ANT <- ggplot(data = linear_rsquared) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$ANT_Imperv_lm, color = 'Imperv')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$ANT_Structure_lm, color = 'Structure')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$ANT_Connect_lm, color = 'Connect')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$ANT_Traffic_lm, color = 'Traffic')) + scale_x_log10() +
labs(title = "Mean ANT ~ R squared of metrics", x ="Focal Distance (meters)", y = "R squared of metrics")
plot(RSquaredPlot_ANT)#graph of R^2 for BIO
RSquaredPlot_Bio <- ggplot(data = linear_rsquared) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$BIO_Imperv_lm, color = 'Imperv')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$BIO_Structure_lm, color = 'Structure')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$BIO_Connect_lm, color = 'Connect')) +
geom_line(aes(x = linear_rsquared$Buffer, y = linear_rsquared$BIO_Traffic_lm, color = 'Traffic')) + scale_x_log10() +
labs(title = "Mean BIO ~ R squared of metrics", x ="Focal Distance (meters)", y = "R squared of metrics")
plot(RSquaredPlot_Bio)